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Mr. E. Bonner was the Program Manager; Drs. W. C. Clever, S. Chakravarthy, 
and V. Shankar served as Pr:lncipal Investigators . 






i 


f 



SUMMARY 


An investigation o£ approximate theoretical techniques for predicting 
aerodynamic characteristics and surface pressures for relatively slender 
vehicles at moderate hypersonic speeds v/as performed. Enphasis was placed on 
approaches that would be responsive to preliminary configuration design level 
of effort. Potential theory was examined in detail to meet 
this objective. 

Numerical pilot codes were developed for relatively simple three 
dimensional geometries to evaluate the capability of the approximate equations 
of motion considered. Results from the computations indicate good agreement 
with higher order solutions and experimental results for a variety of wing, 
body, and wing-body shapes for values of the hypersonic similarity parameter 
M6 approaching one. 
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1. INTRODUCTION 


An examination of the literature for airbreathing hypersonic concepts provides 
an indication of the flexibility and generality required for a prediction 
technique. Typical configuration development variables include wing section, 
incidence, height, dihedral, planform, effectiveness of longitudinal control 
surfaces for trim, effectiveness of en^Dennage for directional stability, and 
propulsion system-airframe interactions. 

State-of-the-art response to these prediction requirements is provided 
by hypersonic impact method*; as well as linearized analysis and design algorithms. 
These approaches can treat significant geometry conplexity with minimum 
response time and cost, witli efficient predicted data coverage in terms of 
Mach number, angle of attack, trim deflection, yaw angle, etc. Shortcomings are 
present, however, in both tlie in 5 )act and linearized methods. For the former, inter- 
ference between surface elements is totally ignored in implementations such as 
classical Newtonian, tangent wedge, and cone theories. Cross -flow interactions and 
stagnation point singularities are also iinplicitly disregarded. In the latter, 
shocks, vorticity, and entropy wakes and layers are excluded. Furthermore, 
superposition of elementary solutions such as those for thickness and angle of 
attack freely used in linear models are, strictly speaking, invalid at hyper- 
sonic speeds. 

A need exists for new aerodynamic prediction techniques to use in the 
optimization of air vehicles designed to travel at hypersonic speeds. Propul- 
sion, structure, and cooling considerations limit the range of cruise Mach 
numbers to be considered to between 4 and 8. One requirement of a new 
aerodynamic prediction technique is that it be more accurate than simple non- 
interfering panel methods. Another specification is that it be more con^juta- 
tionally efficient than currently available explicit finite-difference methods 
so that it can be incorporated into a practical design procedure. The new 
approach should include enough of the physics of the flow to allow realistic 
optimization and should peimit consideration of appropriate interactions 
between conponents of promising hypersonic arrangements, since this has been 
found to be the key to increjasing aerodynamic efficiency at supersonic speeds. 

Less exact non-linear theoretical formulations hold the promise of meeting this 
objective and providing economic design codes which are responsive to 
preliminary vehicle definition efforts. 
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2. LIST OF SmBOLS’^ 


c mean aerodynamic chord 

Cd drag coefficient, D 

CIS 

Cl lift coefficient, L 

qS 

Cm pitching moment coefficient, M 

qSc 

Cp pressixre coefficient, P-Poo 

qco 

D drag 

L lift 

M Mach number or moment 

P static pressure 

q dynamic pressure 

5 reference area 

a angle of attack 

6 flow deflection angle 

SUBSCRIPTS 

<» free stream 


^Additional specialized nomenclature is defined in the theoretical sections 
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3. METHODOLOGY 


Emphasis is placed on approximate theoretical approaches which are capable 
of treating relatively general three dimensional problems but still sufficiently 
sin 5 )le to be responsive to vehicle preliminary design efforts. The basic intent 
of the methodolgy is to produce future improvement in lift-drag ratio of 
hypersonic cruise vehicles. As a result of the strong impact that favorable 
interference has had on supersonic design and the use of such concepts in recent 
advanced hypersonic aircraft studies, candidate analysis should be general enougii 
to systematically treat sucJi problems. Finally, interest in high aerodynamic 
efficiency emphasizes relatively slender configurations at modest angle of attack; 
that is, moderate values of the hypersonic similarity parameter. 

Theoretical effort was recently undertaken (1) to advance hypersonic aerodynamic 
prediction capability at the preliminary design level. A numerical three- 
dimensional second-order refinement of the Prandtl-Glauert model proposed by Van 
Dyke (2) vas developed as a first step up from linear theory. Siich a formulation 
incoiporates nonlinear behavior but retains the isentropic approximimation. This 
approach is known to extend the prediction success for airfoil and cone surface 
pressure to substantially higher values of the hypersonic similarity parameter than 
first-order theory. Typical, three dimensional prediction improvement provided by 
the numerical second-order analysis is presented in • figures 1 and 2 . These results 
indicate sufficient promise to pursue extension of the analysis to treat simple wing- 
body combinations. In addition, the study of second-order potential 
theory indicates that the next level of theoretical richness vis-a-vis a full 
potential equation of motion foimulation should be explored as a means 
of removing edge singularities and ijiproving treatment of characteristic surfaces. 

Hypersonic small disturbance theory was considered in an earlier study (1) 
in recognition of the progressive non- isentropic behavior of the flow as the value 
of the hypersonic similarity parameter increases. Finite difference analysis of this 
approximation (3) indicated that the solution was essentially as complex as that 
for the Euler equation and thus would not be particularly responsive to preliminary 
design level of effort. lids approach was not pursued in the present study on the 
basis of this finding and the previously cited success of potential analysis 
at moderate hypersonic conditions. 
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SECOND ORDER POTENTIAL ANALYSIS 


4.1 APPROACH 

The solution accurate to second order for the flow around wing-body 
cOTibinatiais of aerodynamic interest may be represented by the velocity 
potential, $ , written as the siro of the first and second order velocity 
potentials. 

f Cx, 


Th.e condition of no flow through the boundary requires that 




on the surface of the configuration. The solution , ^ , accurate to first 
order, must satisfy the Prandtl-Glauert equation. 


and the the boundary condition 
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on the surface. For the configurations of interest the aerodynamic surfaces 
are assumed to be thin, and the first order boundary condition is satisfied on 
the mean camber line in accordance with thin wing theory. 
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Such a solution may be obtained by placing source and doublet singularities of 
appropriate strength on or inside the surface of the configuration. 

The method used in this report to obtain first order solutions utilized 
axial singularities and quadrilateral source panels on the body surface in 
conjunction with quadrilateral source and vortex panels to represent the 
(thin) aerodynamic surfaces. The vortex panels, used to represent lift, were 
of constant strength, while the source panels, used to represent thickness. 
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could be made to vary linearly in either the chordwise or spanwise direction. 
The vortex panel and body source panel strengths where obtained by solving a 
set of simultaneous linear equations and thereby satisfying the boundary 
conditions at a set of control points. When the panel singularity strengths 
were obtained the flow properties could be obtained anywhere in the field. 

•Die second order potential, utilizes the first order solution. It must 
satisfy the nonhomogeneous Prandtl-Glauert equation 
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as well as enable the boundary conditions to be satisfied to second order on 
the surface. 


A second order solution using the source volume formulation described in 
reference 1, encounters nintierical difficulties for complex configurations, 
wings with subsonic edges, and supersonic wing-body configurations. The 
primary reason is that source volume strengths require the calculation of 
spatial derivatives, and these cannot be obtained accurately enough from the 
first order solution when using the panel formulation. In addition the 
velocity discontinuities present in supersonic flow are accentuated, since the 
discontinuities introdix:ed by the panel corners and edges do not attenuate 
with distance. Also no reasonable nesh density of spatial source volumes can 
properly account for the large gradients in flow properties near subsonic 
edges. 


Therefore the second order solution was obtained using an approximate 
method which does not require the use of a spatial distribution of sources. 
This solution can be used even when large gradients ace present (e.g. with 
subsonic edges) . The solution is a three dimensional modification of the 
exact second order solution for the pressure coefficient on thin airfoils in 
two dimensions. 






where □ <p * o , 
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-4- f I ^ ^ ] 


and 




'/'(x,y) represents the upper or lower surface- of the airfoil. This solution, 
in effect, represents a local type solution, and all first order velocities 
represented above should be in a coordinate system rotated to make the z-axis 
normal to the local surface (e.g. with nonplanar configurations), but with no 



local angle of incidence. 

Ihe derivation for this method of solution is presented in section 4.3, 
and the derivation for the exact two dimensional solution, from ref 2, is 
presented in section 4.2. Ihe solution uses only the first order flow 
properties on th^ aerodynamic surfaces in conjunction with another first order 
type solution to the homogeneous Prandtl-Glauert equation. Ihis other first 
order type solution is necessary to satisfy the boundary condition to second 
order . 
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A particular solution using axial line sources, accurate to second order, 
was used to help represent the body. 





This particular solution satisfies the second order axisymmetric 
nonhonogeneous Prandtl-Glauert equation to second order. 
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where ^(x,r) is the first order solution to an equivalent axisymmetric body. 
This particular solution was combined with a solution to the homogeneous 
equation. 
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in order to satisfy the boundary conditions. 

The second order accurate axisymmetric axial flow solution was connbined 
with first order cross flow and axial solutions using the body paneling to 



obtain the pressures on the body. The resulting solution is actually a hybrid 
solution combining second order accurate axial velocities with first order 
non-ax isymme trie and cross flow velocities. Although a second order pressure 
coefficient is sufficient for plane type flows, for flows around bodies the 
exact isentropic pressure formula: 
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was found to result in far more accurate pressure predictions when used to 
calculate flow about cones. 

In a sense, the second order requirements on the body are not as 
stringent as on the aerodynamic surfaces since the boundary conditions on the 
body are satisfied on the actual surface, rather than the mean camber line, 
and the cross flow velocities may bfr small compared to the axial velocities. 

In addition the nonlinear expression for Cg in terms of the velocities is used 
on the body. The terms missing from the complete second order solution on the 
body are due to the second order non-ax isymme trie perturbation velocities fran 
the surfaces and from the non-horaogeneous Prandtl-Glauert equation. 


4.2 exact two dimensional ANALYSIS 


Consider flow past a thin airfoil at angle of attack where <x.o 

and <H^are the inclinations of the mean camber line and the free stream with, 
respect to the x-axis. 





TO second order the solution may be written 
where 4 and 4 are solutions to 
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and 
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o 
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with boundary conditions 


^ *Cx) - 


a) , 




To second order the pressure coefficient on the surface is 




To facilitate the solution we define 


4<( 


= J 


'f . 

3t “Z 

'f'» =-0-nl)^S 


and 




<^CXj2) such that 


<P= 




Therefore cD = <D and 
“x* ^zx 


O f = c> 


a <P * o 
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Ihe boundary conditions can be rewritten 


i <x,o) . t - •<, - •<, 
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Ca) 

^ (x,z) may be written as the sum of a particular integral and a solution to 
the homogeneous equation. 

Tex..) / 

Ihe particular integral may be written immediately 


Cx -z) a S- f 


because 
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and at 2 = 0 


<0 M* f ' f 

\ 1^] ^ CVO‘-< h C'f'-'*; -.*]} 


- L ' ’ 1 [ t ^ + <•'*'-% * 1 
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therefore the solution for 4> is complete if 
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and at z ® 0 


U) ^ ‘ IX 


therefore 
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Ihe angle of incidence of tte free stream may be eliminated fron the 
above expression. This means the solution to order two will depend only on 


A » <ai + «< » 

O OO 


Let 


then since 
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and at z * 0 




a> 




O ' ' «K 




U) 
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and therefore 
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^ .. ..,_. r I ^ £!s 2i — if 2, ^ ^ ^ )/*)*• 4> 1 


Since each term in the above expression depends only upon the 

pressure coefficient to sec:ond order is the same whether the free stream or 
the airfoil is placed at an angle of incidence. Replacing y with its 
expression in terms of geometry yields the expression. 
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4.3 EXACT SKEWED WING ANALYSIS , 

A 

Let 4>( S , 1 ) be the first order solution which satisfies 



<^>5(5.0) - 


then Cp to second order may be written; 
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where <^( 5 , ^ ) is a solution to 


O* = o 




and 
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This exact two dimensional solution may be applied to an infinite skewed 
wing in three dimensions because with a rotation of coordinates all flow 
properties depend only on two dimensions. Consider the following 
tr ansf orma tion . 
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The two dimensional solution may be written in three dimensions as follows: 




where 




and 


^ ^ * ax*- ajj* 32 


o <p =* o 




2 ) » O 


However, the exact solution in three dimensions is: 




5»x " -^2 


0"f = O 


n>e only diffoconce betv^en the t«o solutions is , with 


= ! [ I -f ] 

^ (i-n*-) L (l-M*^) J 
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multiplying the second order Cp terms, and the boundary conditions for §• 

The coefficient (1-M« ) multiplying the term is the factor which would be 
present if mass-flux boundary conditions were used, and the second term is due 
to the perturbations of the spatial source distribution. 



^(x,y, 2 ) * 



on the plane 2 * 0 . 

When applied in three dimensions the transformed two dimensional equation 
should give the exact second order pressures vrfienever the flow is locally two 
dimensional and may be approximated by an infinite skewed system, 'flie Mach 
number, M, used in calculating the coefficient is equal to M^cos© (i.e. 
the normal Mach number). Whenever M^jCos© is equal to unity, becomes 
infinite, indicating the second order solution is not valid near M = 1. If 
is used for calculating 5 ^ the solution remains finite except near M^= 1. 

Since is a slowly varying function of M (see figure 4.1.1), except near 
M * 1, (Mg^ ) could be substituted for ^(M,, ) without much change in the 
result. 


4.4 SOLUTION STRUCTURE 


Consider a span station vAiere second order characteristics are to be 
ccmputed, and write all variables in a coordinate system rotated to make the 
2 -axis normal to the local surface. 


z 



Since the first order solution is linear, the first order perturbation 
velocities can be considered to be the composed of four parts. 

1. Contributions from all panels due to an incremental 

change in angle of attack (add load velocities) . 

2. The incremental contribution from all panels due to 

the twist and camber of the aerodynamic surfaces. 

3. The velocity perturbations due to thickness. 

4. The velocity perturbations due to the body in axial flow. 
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Some components have an odd or even part with respect to the local z = 0 
plane. Therefore we can write, 


' V » e<. [V^ t 1- V[v^ii AV^l Ayjl + [ Vj, t AVj.1 

<i>^ z W r«KW^ + r [ 1 -V vJ|^ 

where 1 indicate u^iper or lower surface. The contributions having a i have 
an odd symmetry with respect to the local z = 0 plane. 

The surface displacements i' are composed of three parts. 

1. Displacements due to angle of attack. 

’ 2. Displacements due to twist and camber. 

3. Displacanents due to thickness. 

Therefore we can write. 
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oC f 


)£ i' ± T ir 

X T 


All contributions to the sec:ond order terms are composed of products of first 
order terms. These terms will have either an odd or an even symmetry with 
respect to the local z * 0 plane. The odd synmetry terms in th Cp expression 
contribute to a ACp across the aerodynamic surface while the even symmetry 
terms contribute to an equal Cp on each side. 


q? » -2f 


Cp 




a) 

Cp * 


- 2 [-j- C » - ^ «o»)- I 


and 


AU^^AU.^] + AU.^AU^2 

YU.^+ TU,^] ♦ + -r AU.^3 

± a|o<^U.^AU.^ + X^UtyAUy + T*US.^AU^ ->• 

+ aV [.‘A^AUy + '^yAU,^] •+ Qu.^A\A^+ (A^AIA.^ ] 

> U.^[«K A\A.^+ Y AWjj + AUj^{«< VA.^ •» XU.^ T lA^ ] | 


17 



t 





it la. X r *.■> ■*- 

c< UJ^ •+. uj^ + T [ u3^ 4. 1 

+ 2 oc}^ + 2 otT u:>^ + 2 V t: u?^ 

+ 2 O foC u)^ ■»• T lO 1 

Jb 0< t T 




tf 


<K 


*x 


’■'f'f * * r'-'tiL * L'i'% * %'!'. 1 

* “xx * *xx ^ ’xx "* ‘ix * •‘xx 


with similar terms for 4 • 

Ihe expression is calculated from the second order boundary condition 
for , Ihe odd symmetry second order boundary conditions will be satisfied 
by a hanogeneous solution using source panels, while the even symmetry second 
order boundary conditions will be satisfied using vortex panels and will 
contribute to a net aQ) across the surface. 








where 

‘P.r, = tL<-<] 

= i L ^ ] 


By setting the various coefficients «x., V, or T equal to zero, the 
second order effect of various combinations of angle of attack, camber and 
thickness may be analyzed from the basic solutions. The terms An^and are 
due to the lift due to thickness and are generally small, and in fact are zero 
for planar configurations. The terms AUj^and are the lift induced by the 
body in axial flow and are also in general small. For planar configurations 
or when there is symmetry about the z = 0 plane, the velocity expressions 
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simplify, 
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5. FULL POTENTIAL ANALYSIS 


The full potential equation either in conservative form 
or nonconservative form is frequently used for the treatment 
of transonic flow fields, where the local Mach number in 
general does not exceed '^1.4. However, if the assumptions 
of irrotationality and isentropicity are reasonably valid 
then the full potential equation is expected to yield results 
comparable to Euler equations, even for supersonic/hypersonic 
flow fields. 


The full potential equation can be obtained from the 

continuity of mass relation by assuming the flow to be 
irrotational and isentropic. In the Cartesian system (x,y,z), 
the continuity of mass can be written as 


3 CPU) acpy) 3 (Pw) 

3x 3y 3z 


0 


(5.13 


where p u, pv, and pw represent mass flux in the x, y, and z 
directions, respectively. For an irrotational and isentropic 
flow, the density p, and the velocity components u, v, and w 
are expressed in terms of the velocity potential (j> as 
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u = (|). 


V - (j> 


y 


w * (j> 




(5.2) 


p = 1 - 


1/y-l 


Here, the density and the velocity components are normalized 
with respect to the freestream values. 

Equation (5.1) is said to be in conservation form. The 
nonconservative form of Eq. (5.1) can be obtained in the 
following manner. Differentiating the mass flux terms in 
Eq. (5.1), we get 

p(n^ + Vy + w^) + up^ + vpy, + wp^ * 0 (5.3) 

where 

^ (uu^ * * v™^) 

'’y “ ■ J * ^y * ^™y^ 

p^ = ■ (uu^ * VV^ 4. W^) 

cl 


and a is the local speed of sound expressed as 


a^ = P 


Y-1 




(5.4) 
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Substituting for p , p and p and writing u, v and w in 
terms of 4> in Eq. (5.3), the nonconservative form of Eq. (5.1) 
can be written as 


(a^ - 


i^)(|) + (a^ - <J)^)<f> 


■ 2*xVxy - ‘ ° 


The full potential equation can either be solved in conservation 
form (Eq. (5.1)) or nonconservative form (Eq. (5.5)). In 
general, conservation form is preferred because it allows for 
the correct weak solution (presence of discontinuities like 
shocks and slip surface) to form in the .course of the 
calculation. 

In this report, two different approaches will be described 
to solve the conservative full potential equation. 


5.1 VECTOR APPROACH 

Equation (5.1) contains three unknown variables: u, v 

and w. The density p is expressed in terms of the velocity 

as shown in Eq. (5.2). In the vector approach, the velocities 

/ 

u, V and w will be obtained by solving the following three 
equations: 

(pu)^ + (pv)y + (pw)^ =0 (5.6) 
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0 


(5.7) 


“ u^, « 0 . C5.8) 

Equation (5.6) is the same as Eq. (5.1) and Eqs . (5.7) and 
(5.8) are the irrotationality conditions. These three equations 
can be conveniently represented in a vector form as 

+ Fy + * 0 (5.9) 

where 



pu 


pv 


pw 

E * 

V 

, F » 

-u 

, G = 

0 


w 


0 


-u 


— — 


— 




Equation (5.9) is hyperbolic with respect to the direction 
X if the velocity component u is greater than the speed of 
sound. Assuming u > a, solution to Eq. (5.9) can be obtained 
by marching in the x direction using any stable and consistent 
finite-difference scheme. Equation (5.9) is in Cartesian frame. 
In order to apply the surface tangency condition at the actual 
body location, a body fitted coordinate transformation is 
essential. If we denote the coordinate transformation as 
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a = a(x) ; n = n(x,y,z) ; C * C(x,y,z) (5.10) 


then, the transformed equation corresponding to (5.9) can be 
written in conservation form as 


/s 

+ Fj, + G„ * 0 
a ^ n 


where 


E = 


F * 


G = 


E 

7 


EC + FC + GC 

— r ^ 


EtIx * Fn + Gn^ 

r 


and J is the Jacobian of the transformation 


(5.11) 


(5.12) 


J = . (5.13) 

In this report, Eq. (5.11) is solved using the MacCormack*^ 
predictor-corrector method assuming a to be the time-like 
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hyperbolic direction. At every marching step, the E vector 
is advanced. Before calculating the F and G vector, the 
density p and the velocity components u, v and w must be 
extracted from the E vector. This decoding process is given 
here. 


S.1.1 Decoding Procedure 



(5.14) 


The velocities v and w are nothing but and e^, respectively. 
The component u is obtained from the relation 






(S.15) 


where i denotes the Newton iteration and f(u) is represented as 


f(u) = Au^”^^ + Bu^’^ + C = 0 


where 


A 


2 “ 
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B * - 1) - 1 

C = . 

Once u is computed from Eq. (5.15) the density is given by 

P - ^ . C5.16) 

Usually, two or three iterations are enough for the recursions 
to converge in Eq. (5.15). 

5.1.2 Shock Jinnp Relations for Potential Flow 

In the present approach the outer shocks are treated as 
computational boundaries and fitted as part of the solution 
by satisfying the appropriate shock jump relations across 
them. 

Let the shock be defined as a constant n surface. 

nCx,y,z) » C . (5.17) 

The unit vector normal to the shock is then defined by 
Csee Fig. 5.1) 

N =* (n^i + riyj + (5.18) 
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where i, J and ^ are the Cartesian unit vectors. Let the 
dependent variables u^, v^^, w^, be known on the upstream 
side of the shock. For the discussion in this report, these 
are freestream values. Then, if the shock geometry defined 
by and n_ is known, one should be able to compute 

^2* ^2* ^2 ^2 downstream side of the shock. On 

the other hand, if the density behind the shock, p 2 > is known 
along with two of the direction cosines of the shock normal, 
say riy. and n^, then ^ 2 , '^2 computed. These 

two situations may be termed the direct and inverse shock 
jump problems. 


5 . 1 . 2 . 1 Direct Shock Jump Relations 

hy, r\^, u^, v^, w^, p^ * 1.0 are known. 




= (u^i + V J + w^k) 


(o^^i Pyj P^k) 


IJ ”r 


y/ni + r\t + n 


X 'y 


^s " l^lNl/ai 


p 2 is computed by Newton iteration to satisfy 




. l)p| 


= 0 
s 
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^IN * 

As an example, is assumed unknown 
A . u| - Ujj, 

B - u^Cviny * 
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(5.20) 


c = (vj - 


(w| 


- 




TI, = (- B ± /b^ - AC) /A 

<A> 


(The particular root to be calculated must be specified) 



The vector approach which solves Eq. (5.11) has both 
advantages and disadvantages. The main advantage is that 
it is casted in a convenient conservative form which is 
readily adaptable to any hyperbolic finite difference scheme, 
especially the MacCormack method. The disadvantages are 
that one must store p, u, v and w at each grid point, which 
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means enormous storage requirements for a three-dimensional 
problem (Euler solvers require even more storage than the 
current full potential vector approach due to the addition 
of the more variable pressure) , and the decoding process which 
involves solving a polynomial is very time consuming. 

In order to alleviate the above mentioned disadvantages 
of the vector approach, an alternate method for solving the 
conservative form of the full potential equation is proposed 
in the next section, 

S.2 SCALAR APPROACH 

The method to be described in this section is termed 
scalar approach because the objective of the method is to 
solve for the scalar velocity potential (j> rather than for 
u, V, w and p as in the vector approach. Such a method would 
eliminate the decoding process described in Section 5.1.1 
and also substantially reduce the storage requirement because 
only one variable, (|), need be stored at each grid point. 

Both these factors should contribute toward increasing the 
computational efficiency and speed. In what follows, an 
implicit marching scheme utilizing a mathematically sound 
density linearization method which solves for the scalar 
velocity potential (|> is described. 

Referring to Fig. 5.2, let us assume that the information 
is given at all points in planes i, i-1, and i-2. The objective 
is to compute the (J> information in plane i+1 by solving Eq. (5.1) 
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1-2 PLANE i-1 PLANE i PLANE i + 1 PLANE 


Figure 5.2. Inplicit Conputational Molecule 
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The implicit finite differenced form of Eq. (5.1) can be 
written as 






l.j.k“l,J,k1 

L 


(5.21) 


0 . 


In Eq, (5.21), the flow quantities like density p and Velocity components 

Uf Vf and w appear at the i+1 level, \»hlch makes it difficult to solve 
for (t> at i+1 without introducing some kind of a linearization. Since 
p is a complicated function of (j), the logical variable to linearize is 
density. 

Consider the first bracketed term in Eq. (5.21) which represents the mass 
flux balance in. the x-direction« To write this in terms of values, 

let us linearize In terms of and p^ in terms of To maintain 

the conservation law form both and p^ have to be linearized. Thus, 


u 


i+l,j,k''i+l,j,k 


Q +1 Ax ' 

Pi,j,k ^ \3x n 


4* •••> u 


i,j,k 


i+1, j,k (5.22) 


where 


|£ » -p^“'^M^(uu^ + + ww^) . 


(5.23) 
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Substituting Eq. (5.23) into Eq. (5.22) and rearranging the teims results 
in the following. For convenience, the subscripts j and k are suppressed. 


^i+l\+l * [^i “ “i+1 


+ + u^v| + u^w|) 




where 


(5.24) 




i+l.j.k 


/wih\ 

I ^ /j.i 


backward differenced 


(5.25) 


w 


i+l.j.k 


*l+l.J+l,k - ♦n-l.j-l.k 

2Ay 


2Az 


central differenced 


Similarly, expanding Pi,j.k about C®ven through Pi,j,k is 
known, to maintain conservative differencing it has to be expanded 
about Pi-i,j^k)» write (suppressing j and k indices) 
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= [pi-1 - 




( 5 . 26 ) 






Combining Eqs. (5.24) and (5.26), one can write the first term in 
Eq. (5.1) as 


3pu . 
3x ~ 


X 


(P. 




i+1 


.2-Y 






-Y 9 

« 1 i i+1 




(5.27) 


In Eq. (5.27), only the velocities u, v, and w appear at the 
unknown i+1 level. The density is linearized in such a way (using 
Eqs. (5.22) and (5.23) that it appears only at the known i and i-1 
levels. Substituting for u, v, and w in terms of the velocity 
potential from Eq. (5.25), and defining a differential velocity 
potential A4> such that A4> * (|). . . , - (j>. . , Eq. (5.27) can be 

l+i,J,K i,j ,K 

finally written into delta form as 
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3(Pu) 

3x 




n2-y 2 

S^i “-"i”! 


3z 


A(j) - 


JL 

Ax 


-2-Y 2„3( 
Pi ^oo“i 


(5.28 


There are some interesting physics featured in Eq,. (5.28) which will 
be discussed later. Nov, consider the second term in Eq.'(5.1) whose 
finite differenced form is given by the second bracketed term in Eq.i (,5.21). 
Here again, both the density and velocity v appear at the i+l level and 
some kind of a density linearization is essential to solve for (() at 
the i+1 level. A simple linearization of the type 


^i+l,j+»l.k “ ‘^i,;rt,k 
^i+l,j-Ji,k “ *^i,j-Js,k 


(5.29) 


is sufficient for the y-derivative treatment in the Cartesian framework 
since the velocity component v is usually subsonic for a predominantly 
supersonic flow in the x-direction. If the v component does go from 
subsonic to supersonic or vice versa, then the density linearization 
given by Eq. (5.29) will have to be modified to a slighy.y complicated 
form, which will not be discussed here. Using Eq. (5.29), the second term 
in Eq. (5.1) or Eq. (5.21) can be written in the delta form as 
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^ ■ W ‘’i+1 ¥ *i+l “ F ‘’l W ^'^1+1 ■ '"l * 1 ’ 


3 . 3 


^ ^ A({. + ^ (p^Vi) . 


(5.30) 


Similarly, the z-derivative term in Eq. (5.1) can be linearized and 
written as 


lE“ i 3 p 3 w + 3 (p „ ) 
3z 3z *^i 3z ^ Jz '•‘^i i*! 


(5.31) 


Combining Eqs. (5.28), (5.30), and (5.31), we have the final density 
linearized foim of the full potential equation. 


Ax= 




+ •>;“ p Ad) + p. A(b “ R 
3y ^i 3y ^ 3z '^i 3z ^ i 


■(5.32) 


where 




i Ax I i i <» i 
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In Eq (5.32), both the y and z derivative terms appear in an implicit 
manner which makes It difficult to solve without some approximations to 
the left hand side. Factoring the left hand side into (x,y) terms and 
(x,z) terms for computational implementation, we get 


^ «*'^i'*^i 3z 3z ^i 3z ^ ^ ^7 3y 3y y ^ i 


■V“ 

A 


-/ V. 


■V 


(5.33) 


In the above eqtiation A and B are operators where the operator A gives 
rise to an implicit molecule in the x-z plane, while B results in an 
implicit molecule in the x-y plane. The solution procedure for Eq. (5.33) 
is twofold. Define A4)* • BA(|) and then perform the following: 


AA4>* - 
BA(j) - A4>* 


(5.34) 


♦ 


i+l» j »k 




The implicit operator A when written down results in a trldiagonal form 
as shown below. 
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( 5 . 35 ) 



For the solution procedure to converge, it is necessary and sufficient 
that the above tridiagonal matrix possess diagonal dominance such that 


bj > |aj + jc| . 


( 5 . 36 ) 


The diagonal term b looks like 


\ 


( 5 . 37 ) 


In order to ensure diagonal dominance the term B/Ax^ has to be negative. 
Let us. go back and review the term S. with the physics of the problem 
in mind. 


2-Y«2„2 

<=1 


( 5 . 38 ) 


Also, the density is related to the speed of sound and Mach number through 
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(5.39) 


Pi »i«« 


where is the local speed of sound. 

Substituting Eq. (5.39) into Eq. (5.38) we get 



(5.40) 


If we assume the flow is supersonic in the marching direction x, then 
it implies > a^, which means. 



(5.41) 

/ 

or 

B < 0 . 

This diagonal dominance property was achieved only through the particular 
form of the density linearization assumed in Eqs. (5.22) and (5.23). 


If we had assumed at that point that p. . . j i, = Pj ^ ^ without the 

-r^ I , then the final expression for B would have resulted 
in B * P^ which would always be positive, and diagonal dominance is not 
ensured. Thus, the form of density linearization is very critical for 
convergence of the solution procedure. 
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5.2.1 Arbitrary Coordinate System 

The entire analysis described in Section 5.2 assumed a 
Cartesian framework where it is convenient to explain the 
density linearization concept and also simple to construct 
a computer program to test the concepts. However, for the 
treatment of wing-body configurations, a body-fitted 
coordinate transformation is essential for application of 
body boundary conditions. In this section, a brief introduction 
is made for the inclusion of a coordinate transformation into 
Eq. (5.1) and the associated density linearization procedure. 

Considering the general transformation given by 
Eq. (5.10), one can write Eq. (5.1) in the form 



(5.42) 


where 


U * ua„ 

Jv 


V = UC 


X 




. vKy * . 
+ vriy + wn^ * + 
1 - ^ 


A«4)j. + A-(^ 

hh * 

+ - 1 ) 




> (5.43) 


42 


A, - 




+ a- 


An *crE + G E G E 
2 x^x y^y z^z 




A-®* an +an +an 
3 X 'x y 'y z 'z 


A4 * A2 


A- » 

5 ^x ^y ^z 


^6 “ + Syny * 


A? - Aj 


^8 ■ h 


A„ = 


n^ + + «2 

'x 'y 'z 


> (5.43) 


In Eq. (5.42), the direction a can be considered as the 
inarching direction (time-like hyperbolic direction), provided 
the projection of the velocity in the direction normal to the 
constant d line/surface is supersonic. This is illustrated 
in Fig. 5.3 and the verification of the above statement will 
become clear once the density linearization is introduced 
into the' first term of Eq. (5.42). 
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Consider the first term in Eq. (5.42). The objective 
now is to linearize the density term p in such a way we can 
not only show the requirement for a to be the marching 
direction but also be able to solve for the velocity poten- 
tial ((>. Let p be expanded about Pq where Pq is some neigh- 
boring known state. 

P “ PO * * •” 


where 


If = - if K ^ 


Using Eq. (5.44) we can write 


/pU\ 

\ZJ 

”3a 


-L 

da 


4jVa ^ Vc " Vn * ^0 


u 

J 


(5.45) 


where 


'o-"o(‘.), -"oWo 


(^) ^ ^(‘=0 j ) _ 


a /i 


aa aa 


^ \7 * " 0 ^ 0*5 * ™ 0 « 0 ^ " 

CS.46) 


4'5 



Substituting for U in terms of A^, A 2 and A^, and rewriting 
Eq. (5.46) , we get 


30* 


a 

30' 


AiPq 


^0 


30 




—pl 


a 


0 


+ ■ 


Vo 


a 


0 / 




f u 
- 0 0 

0 ‘ „2 

^0 


(5.47) 


where 


3 ( ) _ ^ j »k ” ^ ^i» j ,k . backward 

3o ” Ao * differenced 


3( ) _ ^ ^i»l>j+l,k ~ ^ ^i^^l,j’l,k 


3C 


1A5 


3( ) ^ ^ ^i-»l,j,k+l ~ ^ ^i-«-l,j,k~l 
3n 2 Ati 


y (5.48) 


central 
differ- 
enced 


To ensure diagonal dominance for stability, the term 
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appearing in Eq. (5.47) has to be negative. This implies 
(i£ backward differencing is used in a) 




or 




; for stability . 


(5.49) 


The quantity is nothing but the projection of the 

velocity vector q in the direction normal to the constant a 
surface as shown in Fig. 5.3. Thus, it is verified that for 
a to be the marching direction, has to be supersonic. 

The Cartesian framework described in Section 5.2 (Eq. (5.27)) 
is contained in Eq. (5.47) as a special case (set o»x; C*y; 
n=z; A^=l; - Aj - 0; J«l; Uo=Uo; . 

The formulation for the scalar approach in the transformed 
coordinate system a,C,n is at present in the conceptual stage 
as far as the treatment of and terms in Eq. (5.42). 

Results from the scalar approach to be presented in this report are 
based on the Cartesian system (Eq. (5.1)). 
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6.0 RESULTS 


6.1 SECOND ORDER COMPARISONS 

Coirparisons were performed to determine the success of the proposed second 
order analysis technique in prediction of surface pressures, forces, and moments. 
The following sections present comparisons between first and second order 
analysis, experimental data, and higher order numerical calculatiions . 

6.1.1 Wing Alone 


Figures 6.1.1 through 6.1.3 present first and second order results for wing 
alone cases. The predictions are based on using perturbation velocities obtained by 
placing the free stream at angle of attack. The comparisons indicate that even for 
wings with subsonic edges and in cases where the flow was far from two dimensional, 
the predictions were quite good. The integrated lift and moment calculation shown 
in figure 6.1.3 demonstrates a substantial improvement in the aerodynamic center 
prediction for a wing with a siibsonic edge. Reference 1 presents similar results 
for a wing with a si; 5 )ersonic leading edge based on a source volume method. The 
present analysis gives virtually the same answer. 

6.1.2 Bodies of Revolution 


The pressures and velocities on bodies of revolution were con^juted by combining 
a first or second order axial solution with a first order cross flow solution. Cone 
thickness and angle of attack con^jarisons between hybrid second order analysis and the 
conplete Euler equations, figures 6.1.4 through 6.1.6, demonstrates that a second 
order analysis can yield considerable inprovement over first order results even for 
values of M6 near 1. The-resiiLts for a non-lifting parabolic body of revolution, 
presented in figure 6.1.7, show that the difference between a first and second 
order analysis is greatest in regions of highest flow deflection angle near the 
nose and tail as would be expected . 

Figures 6.1.4 and 6.1,5, demonstrate that the second order velocities computed 
for axial flow past cones conpares very favorably with the exact solutions . The 
second order solution appears to be accurate to values of M6 approaching 1. . 

Figure 6.1.6 show how the isentropic pressure coefficient formula, when combined 
with the accurate velocity predictions resiilt in good prediction of the pressure 
coef f icients - for lifting cones. 

Results for M == 6 flow past the body shown in figure 6.1.8 are presented 
in figure 6.1.9. For analysis purposes the nose region of the body has to be 
modified to avoid slopes which exceed the Mach angle. A conical forebody extension 
having a slope slightly less than the Mach angle was used for this purpose. In the 
nose region, where the values of M(S were large, differences between the analysis and 
the experimental data are noticeable. The greater slopes present on some bodies 
can restrict regions which may be usefully analyzed. 



6.1.3 Delta Wing-Body 


Figure 6.1.10 present results for the pressure coefficient calculation for 
flow past the wing-body coihbijiation shown in figure 6.1.18:. For most regions on the 
wing, a comparison with the €>xperimental data is reasonably good. 
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Figure 6.1.4. Supersonic Thickness Velocities for 10 Degree Cone 














































































































c) a variable, ,0 = 90 ‘ 


Figure 6.1.9. Conpleted 


































6.2 FULL POTENTIAL COMPARISONS 

In this section, results will be presented for both the 
full potential vector approach (Section 5.1) and the full 
potential scalar approach (Section 5.2). To verify the 
validity of the full potential theory, results are presented 
for various unit problems such as supersonic flow over 
two-dimensional wedge, supersonic leading edge delta wing, 
subsonic leading edge delta wing, and a simple conical 
wing-body combination. 

6.2.1 Conical Wedge Flow 

Figure 6.2.1 shows the schematic diagram of supersonic 
flow past a wedge with an. attached shock wave. To implement 
the vector approach, a conical transformation of the independent 
variables was chosen to align the shock with a constant 
coordinate line. 




This transformation coupled with the use of the strong 
conservation law form of the vector equation (5.11), and the 
conservative differencing provided by the MacCormack scheme, 
results in fairly weak oscillations across the shock. 

To start the solution procedure, freestream values are 
specified at all grid points at an initial x location 
(x * 1.0). The solution is then marched along the x direction 
until an asymptotic steady state is reached. Care is taken 
to make sure that the grid points chosen enclose the shock. 
Freestream values are maintained at the outer boundary while 
surface tangency is imposed on the wedge boundary. 

A comparison of the solutions obtained by using the Euler 
equation, the vector approach full potential equation (Eq. (5.11)), 
and the Cartesian, scalar approach (Eq. (5.33)), is shown in 
Fig. 6.2.2. The results from the Euler and full potential 
vector approach show noticeable oscillations across the shock 
wave while the scalar approach seems to provide smooth monotone 
results. For this simple supersonic wedge flow example, the 
scalar approach demonstrated at least a factor of 10 speed-up 
in the computational time over the Euler solver as well as 
the vector approach. Also shown in Fig. 6.2.2 are the 
theoretical linear and second order solution. 
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6 . 2.2 


Supersonic Leading Edge Delta Wing 
For a flat plate delta wing with a supersonic leading 
edge, the solution on the windward and leeward sides are 
independent of each other. We will study the windward or 
compression side in this report for this will test the 
numerical method’s ability to capture three dimensional shocks 
and also check out the equations developed to fit the leading 
edge shock (see Fig. 6.2.3 for a schematic view of the flow 
field). For both shock capturing and shock fitting calculations, 
an initial solution is specified on an initial data plane at 
X » 0.1. The solution is then advanced downstream until it 
converges to a conically similar steady state. 

6-. 2 . 2 . 1 Shock Capturing 

The coordinate transformations used in the calculations 
where the shock is captured are 

0 * X 

C * z/x (6.2.2) 

p = y/ (x - z tan A) 

The constant coordinate lines in the physical and computational 
domains are sketched in Fig. 6.2.4. Freestream values are 
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LEADING EDGE 



Figure 6.2.3. Schematic Diagram o£ Delta Wing Flow Field 
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specified at all grid points in the initial data plane. For 
this coordinate system it is advantageous to specify zero 
derivatives in the C direction at the boundary BC where the 
flow is two-dimensional, 

6 • 2 . 2 . 2 Shock Fitting 

The coordinate transformations used in the calculations 
where the shock is fit are given by 

CT * X 

C - z/x (6.2.3) 

„ = ^ 

^shock " ^body 

The corresponding constant coordinate lines are sketched in 
Fig. 6.2.5. The exact two-dimensional solution is imposed on 
the boundary BC at every marching step. The shock jump 
relations are imposed at the shock boundary CD. To start 
the calculations, an initial shock shape is chosen which 
determines the solution at all grid points on CD. At every 
spanwise grid point the initial values of the dependent 
variables are assumed to be those at the corresponding shock 
grid point. The shock slope at C is fixed at the exact value 
for all steps. 
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For the vector approach, both the shock capturing and 
shock fitting schemes were implemented while the scalar 
approach results presented here are based on shock capturing 
procedure. In Fig. 6.2.6, pressure coefficient results from 
the full potential vector and scalar approach are compared 
with results obtained with a shock capturing Euler solver and 
second order potential theory for * 4.0, A * 50® and 
a * -5®, -10® and -15®. Although the full potential scalar 
approach and vector approach predict identical results, the 
scalar approach is an order of magnitude faster than the 
Euler or the full potential vector approach. For low angles 
of attack Ca = -5®), the full potential solution is 
indistinguishable from the Euler solution. 

From comparing the second order theory solution with the 
FPE and the Euler solutions, the inability of the second order 
theory implementation to shift the cross flow sonic line 
Cjunction of the two-dimensional and three-dimensional regions) 
becomes clear. The FPE is superior in this respect, as far as producing 
results comparable to Euler solvers. 

6.2.3 Subsonic Leading Edge Delta Wing 

One of the advantages of full potential theory 
over the linear or second order theory is that it does not 
create any leading edge singularities, and thus can easily 
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Figure 6.2.6. Conparison o£ Pressure Coefficient on Supersonic Delta Wing 
at Various Angles of Attack, M» = 4.0, A = 40° 





handle subsonic leading edges. The second order finite difference 
results for the subsonic leading edge cases are shown in 
Ref. 1 which indicate pronounced oscillations in pressure 
in the neighborhood of the leading edge because of the 
singularity inherent in the theory. This behaviour 
is avoided by using full potential theory as illustrated 
in Figs. 6.2.7 and 6.2.8. Figure 6.2.7 shows the results 
obtained using the full potential vector approach. All four 
cases iX « 20®4», 28*50’, 38*39’ and 43*38’) are presented. 
There are no oscillations in the surface pressure coefficient. 
Some breaks in smoothness are evident near the leading edge 
on the wing surface. Tests reveal that this is due to the 
discontinuity in the grid used (analytic metrics were also 
used for these calculations and so the mesh was not implicitly 
smoothed out by the numerical evaluation of metrics) . For 
the low freestream Mach number of 2 considered, tests showed 
that the solutions to the FPE and Euler equations are 
indistinguishable because of the negligible entropy generation. 
However, this would be true even for higher Mach numbers for 
slender delta wings with subsonic leading edge because the 
leading edge compression system is relatively weak for these 
cases except possibly when the leading edge is very close to 
the Mach cone associated with the freestream Mach number. 
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Figure 6.2.7. Pressure Coefficient on Series of Delta Wings with Subsonic 
Leading Edges using the Full Potential Vector Approach 








The results from the full potential scalar approach are 
compared in Fig. 6.2.8 with the vector approach results. The 
scalar approach results are much superior and also required 
much less time to compute over the vector approach. 

6.2.4 Cosine Wing-Body Combination 

The wing -body configuration consisted of a flat plate 
delta wing with a conical fuselage represented by the formula 

^ ° < I < “-2098 . 

The shock capturing versions of the finite difference 
methods using the Euler equation and the full potential vector 
and scalar approaches were used to compute the pressure 
coefficient on the surface of the wing-body combination 
and the results are shown in Fig, 6.2.9, It is evident that 
for the low angle of incidence considered here the full 
potential and Euler solutions are in good agreement. The 
scalar approach program utilizes a linearization of the actual 
body boundary conditions because the code at present is in 
terms of a Cartesian system in which the body is not in 
general aligned with grid points. 

The results are in qualitative agreement with the results 

5 

obtained by Kutler for a similar geometry. The salient 
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features of the pressure distribution on wing-body combinations 
similar to the ones considered here are as follows (1) A 
two-dimensional region between the supersonic leading edge 
and cross-flow shock (a part of which is shown as AB in 
Fig. 6.2.5). (2) A pressure rise across the cross-flow shock 

denoted by BC. (3) Compression up to the junction of wing 
and body denoted by CD. (4) An expansion followed by compression 
on the fuselage denoted by DE. 

Here again the scalar approach demonstrated a significant 
Cat least a factor of 10) improvement in computational speed 
over the Euler solver or the full potential vector approach. 
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7. CONCLUSIONS 


Based on the theoretical development and comparison with higher order 
results/experimental measurements, the following conclusions are made. 

1. Inproved prediction of si 5 }ersonic/hypersonic aerodynamic characteristics 
and surface pressures for a single wing, body and wing-body shapes has 
been demonstrated for non-linear potential analysis. 

2. Second order theory provides a systematic means of extending linear 
analysis to values of the similarity parameter M6 approaching one. 

Fifty second CPU solution time/Mach number is typical for a wing- 
body problem. 

3. Full potential analysis successfully eliminates subsonic edge 
singularities and linear characteristic approximations of second 
order theory. The scalar formulation is an order of magnitude 
faster than the vector approach and Euler solvers. 

4. Potential theory provides an advanced aerodynamic prediction technique 
that is responsive to the preliminary design problem at moderate 
hypersonic conditions . 

5. Further effort is required to extend the analyses to more general 
geometries and develop configuration design codes. 
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